libraries
library(tidyverse)
library(brms)
library(bayesplot)
library(here)
library(glue)
library(scales)
theme_set(theme_light())
source(glue("{params$common_dir_str}/brms_model.R"))
source(glue("{params$model_dir_str}/model_prior.R"))
data
obs_only <-
read_csv(glue("{params$model_dir_str}/data/stimulation_obvs.csv")) %>%
mutate(subj = as_factor(subj),
obs_degree = error,
error = obs_degree * (pi/180))
## Parsed with column specification:
## cols(
## subj = col_double(),
## subj_index = col_double(),
## stimulation = col_double(),
## error = col_double()
## )
full model
print(bprior_full)
## prior class coef group resp dpar nlpar bound
## 1 normal(3.8, 0.4) b intercept circSD
## 2 normal(0, 0.4) b stimulation circSD
## 3 normal(0, 0.25) sd Intercept subj circSD
## 4 normal(0, 0.25) sd stimulation subj circSD
## 5 normal(0, 1.5) b intercept theta
## 6 normal(0, 1) b stimulation theta
## 7 normal(0, 0.5) sd Intercept subj theta
## 8 normal(0, 0.5) sd stimulation subj theta
iter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains
fit_full <- brm(bform_full, obs_only, prior = bprior_full,
family = vm_uniform_mix, stanvars = stanvars,
warmup = warmup, iter = iter, cores = cores, chains = chains,
control = list(adapt_delta = 0.99), inits = 0,
file = glue("{params$save_dir_str}/fit_full"))
## Compiling the C++ model
## Start sampling
fit check
divergences
model_fit <- fit_full
#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_div = sum(Value))
## n_div
## 1 0
rhat
mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neff ratio
mcmc_neff(neff_rat) + yaxis_text(hjust = 1)
trace plots
mcmc_trace(as.array(model_fit$fit))
reduced model: no pre to post change on circSD parameter
print(bprior_DcircSD_null)
## prior class coef group resp dpar nlpar bound
## 1 normal(3.8, 0.4) b intercept circSD
## 2 normal(0, 0.25) sd Intercept subj circSD
## 3 normal(0, 0.25) sd stimulation subj circSD
## 4 normal(0, 1.5) b intercept theta
## 5 normal(0, 1) b stimulation theta
## 6 normal(0, 0.5) sd Intercept subj theta
## 7 normal(0, 0.5) sd stimulation subj theta
fit_DcircSD_null <- brm(bform_DcircSD_null, obs_only, prior = bprior_DcircSD_null,
family = vm_uniform_mix, stanvars = stanvars,
warmup = warmup, iter = iter, cores = cores, chains = chains,
control = list(adapt_delta = 0.99), inits = 0,
file = glue("{params$save_dir_str}/fit_DcircSD_null"))
## Compiling the C++ model
## Start sampling
fit check
divergences
model_fit <- fit_DcircSD_null
#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_div = sum(Value))
## n_div
## 1 0
rhat
mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neff ratio
mcmc_neff(neff_rat) + yaxis_text(hjust = 1)
trace plots
mcmc_trace(as.array(model_fit$fit))
reduced model: no pre to post change on pMem parameter
print(bprior_DpMem_null)
## prior class coef group resp dpar nlpar bound
## 1 normal(3.8, 0.4) b intercept circSD
## 2 normal(0, 0.4) b stimulation circSD
## 3 normal(0, 0.25) sd Intercept subj circSD
## 4 normal(0, 0.25) sd stimulation subj circSD
## 5 normal(0, 1.5) b intercept theta
## 6 normal(0, 0.5) sd Intercept subj theta
## 7 normal(0, 0.5) sd stimulation subj theta
fit_DpMem_null <- brm(bform_DpMem_null, obs_only, prior = bprior_DpMem_null,
family = vm_uniform_mix, stanvars = stanvars,
warmup = warmup, iter = iter, cores = cores, chains = chains,
control = list(adapt_delta = 0.99), inits = 0,
file = glue("{params$save_dir_str}/fit_DpMem_null"))
## Compiling the C++ model
## Start sampling
fit check
divergences
model_fit <- fit_DpMem_null
#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_div = sum(Value))
## n_div
## 1 0
rhat
mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neff ratio
mcmc_neff(neff_rat) + yaxis_text(hjust = 1)
trace plots
mcmc_trace(as.array(model_fit$fit))
reduced model: no pre to post change on pMem or circSD parameter
bprior_DcircSD_DpMem_null
## prior class coef group resp dpar nlpar bound
## 1 normal(3.8, 0.4) b intercept circSD
## 2 normal(0, 0.25) sd Intercept subj circSD
## 3 normal(0, 0.25) sd stimulation subj circSD
## 4 normal(0, 1.5) b intercept theta
## 5 normal(0, 0.5) sd Intercept subj theta
## 6 normal(0, 0.5) sd stimulation subj theta
fit_DcircSD_DpMem_null <-
brm(bform_DcircSD_DpMem_null, obs_only, prior = bprior_DcircSD_DpMem_null,
family = vm_uniform_mix, stanvars = stanvars,
warmup = warmup, iter = iter, cores = cores, chains = chains,
control = list(adapt_delta = 0.99), inits = 0,
file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))
## Compiling the C++ model
## Start sampling
fit check
divergences
model_fit <- fit_DcircSD_DpMem_null
#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)
np %>%
filter(Parameter == "divergent__") %>%
summarise(n_div = sum(Value))
## n_div
## 1 0
rhat
mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
neff ratio
mcmc_neff(neff_rat) + yaxis_text(hjust = 1)
trace plots
mcmc_trace(as.array(model_fit$fit))
loo compare
#need to do this only once?
expose_functions(fit_full, vectorize = TRUE)
fit_full <- add_loo(fit_full, file = glue("{params$save_dir_str}/fit_full"))
fit_DcircSD_null <- add_loo(fit_DcircSD_null, file = glue("{params$save_dir_str}/fit_DcircSD_null"))
fit_DpMem_null <- add_loo(fit_DpMem_null, file = glue("{params$save_dir_str}/fit_DpMem_null"))
fit_DcircSD_DpMem_null <- add_loo(fit_DcircSD_DpMem_null, file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))
print("fit_full")
## [1] "fit_full"
fit_full$loo
##
## Computed from 8000 by 504 log-likelihood matrix
##
## Estimate SE
## elpd_loo -805.2 15.9
## p_loo 6.7 0.2
## looic 1610.4 31.7
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
print("fit_DcircSD_null")
## [1] "fit_DcircSD_null"
fit_DcircSD_null$loo
##
## Computed from 8000 by 504 log-likelihood matrix
##
## Estimate SE
## elpd_loo -805.5 15.8
## p_loo 6.3 0.1
## looic 1610.9 31.6
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
print("fit_DpMem_null")
## [1] "fit_DpMem_null"
fit_DpMem_null$loo
##
## Computed from 8000 by 504 log-likelihood matrix
##
## Estimate SE
## elpd_loo -804.5 15.9
## p_loo 6.1 0.1
## looic 1608.9 31.8
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
print("fit_DcircSD_DpMem_null")
## [1] "fit_DcircSD_DpMem_null"
fit_DcircSD_DpMem_null$loo
##
## Computed from 8000 by 504 log-likelihood matrix
##
## Estimate SE
## elpd_loo -804.9 15.8
## p_loo 6.0 0.1
## looic 1609.9 31.6
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
loo_compare(fit_full, fit_DcircSD_null, fit_DpMem_null, fit_DcircSD_DpMem_null, criterion = "loo")
## elpd_diff se_diff
## fit_DpMem_null 0.0 0.0
## fit_DcircSD_DpMem_null -0.5 0.8
## fit_full -0.8 0.4
## fit_DcircSD_null -1.0 1.0